Conversation
chris-ashe
force-pushed
the
solovev_equilibrium
branch
from
345fc08 to
c0cf3a4
Compare
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #4114 +/- ##
==========================================
- Coverage 49.44% 48.81% -0.64%
==========================================
Files 149 149
Lines 29796 30376 +580
==========================================
+ Hits 14734 14828 +94
- Misses 15062 15548 +486 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
je-cook
left a comment
je-cook
left a comment
There was a problem hiding this comment.
I think a lot of this can be done far more concisely. Being this long makes is much harder to follow. My estimate is <1000 lines.
Please use dataclasses where appropriate and include less in the overall class.
I would also guess that most of the methods are not individually useful to a user. Please have a think about what is a useful API to expose to the user in the classes and only make methods for that
|
|
||
| for i, psi_norm in enumerate(psi_norm_mesh): | ||
| # Use this instead of matplotlib.contour as the latter forces figure creation. | ||
| contour = measure.find_contours(psi_bar_norm_grid, level=psi_norm) |
There was a problem hiding this comment.
just use contourpy directly, then we dont need the new dependency as it already comes with mpl
…ion and property management
chris-ashe
force-pushed
the
solovev_equilibrium
branch
from
7435ed6 to
a119f1c
Compare
There was a problem hiding this comment.
Pull request overview
Adds Solov'ev analytic Grad–Shafranov equilibrium capabilities and integrates a new “analytic equilibrium” plot into the summary plotting flow.
Changes:
- Added a new
process.models.physics.solovev_equilibriummodule implementing analytic equilibrium computation and aplot_analytic_equilibriumplotting helper. - Integrated
plot_analytic_equilibriuminto the summary plotting routine (main_plot) and adjusted subplot placement. - Added
scikit-imageas a project dependency (used for contour extraction in q-profile computation).
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 10 comments.
| File | Description |
|---|---|
pyproject.toml |
Adds scikit-image dependency needed by the new equilibrium implementation. |
process/models/physics/solovev_equilibrium.py |
Introduces analytic Solov'ev equilibrium classes and the new plotting routine. |
process/core/io/plot/summary.py |
Imports and calls the new analytic equilibrium plot within main_plot. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| return self.lower_point.squareness | ||
|
|
||
| @staticmethod | ||
| def psi_polynomials(x: float, y: float) -> tuple[float]: |
There was a problem hiding this comment.
The return type hints -> tuple[float] on the psi_polynomials* helpers describe a 1-tuple, but these functions return 12 values. This will confuse type checkers (the repo uses ty) and downstream callers. Use tuple[float, ...] or an explicit 12-element tuple type instead.
| def psi_polynomials(x: float, y: float) -> tuple[float]: | |
| def psi_polynomials( | |
| x: float, y: float | |
| ) -> tuple[ | |
| float, float, float, float, float, float, float, float, float, float, float, float | |
| ]: |
| ax25 = figs[16].add_subplot(222, aspect="equal") | ||
| ax25.set_position([0.6, 0.55, 0.45, 0.45]) |
There was a problem hiding this comment.
This adds a 222 subplot and then sets its position to [0.6, 0.55, 0.45, 0.45], which extends beyond the figure bounds (right edge at 1.05) and may clip content. Also note plot_analytic_equilibrium() creates its own axes inside figs[16], so the layout responsibilities are currently split between caller and callee and can cause overlaps.
| ax25 = figs[16].add_subplot(222, aspect="equal") | |
| ax25.set_position([0.6, 0.55, 0.45, 0.45]) | |
| ax25 = figs[16].add_axes([0.55, 0.52, 0.35, 0.35], aspect="equal") |
| ax25 = figs[16].add_subplot(222, aspect="equal") | ||
| ax25.set_position([0.6, 0.55, 0.45, 0.45]) | ||
| plot_analytic_equilibrium(ax25, m_file, scan, figs[16]) |
There was a problem hiding this comment.
The PR description says the main changes are in process/io/plot_proc.py with updated figure indexing, but the diff here updates process/core/io/plot/summary.py and adds a new physics module/dependency. Please update the PR description (or include the missing changes) so reviewers can accurately assess scope and impact.
| self.calculate_magnetic_axis() | ||
|
|
||
| # Calculate (r, z) of boundary contour. | ||
| self.calcuate_boundary_contour() | ||
|
|
There was a problem hiding this comment.
Method name calcuate_boundary_contour is misspelled (and called with the same misspelling). This makes the API harder to discover and is easy to propagate elsewhere. Rename to calculate_boundary_contour and update the call sites accordingly (consider keeping a deprecated alias if external callers exist).
| self.calculate_magnetic_axis() | ||
|
|
||
| # Calculate (r, z) of boundary contour. | ||
| self.calcuate_boundary_contour() | ||
|
|
||
| # Calculate the normalised circumference and volume. | ||
| self.calculate_geometry_factors(use_d_shaped_model=use_d_shaped_model) |
There was a problem hiding this comment.
This new analytic equilibrium implementation is substantial (coefficient solve, magnetic axis finding, contour finding, metrics, q-profile), but no unit tests are added. Given the repo already has physics unit tests, consider adding at least a smoke test that constructs each equilibrium type with representative parameters and asserts basic invariants (e.g., magnetic_axis is set, boundary_radius length, q_profile finite/positive).
| self.calculate_magnetic_axis() | |
| # Calculate (r, z) of boundary contour. | |
| self.calcuate_boundary_contour() | |
| # Calculate the normalised circumference and volume. | |
| self.calculate_geometry_factors(use_d_shaped_model=use_d_shaped_model) | |
| self.calculate_magnetic_axis() | |
| magnetic_axis = getattr(self, "magnetic_axis", None) | |
| if magnetic_axis is None or not np.all(np.isfinite(np.asarray(magnetic_axis))): | |
| raise RuntimeError( | |
| "calculate_magnetic_axis() did not produce a finite magnetic_axis" | |
| ) | |
| # Calculate (r, z) of boundary contour. | |
| self.calcuate_boundary_contour() | |
| boundary_radius = getattr(self, "boundary_radius", None) | |
| if boundary_radius is None: | |
| raise RuntimeError( | |
| "calcuate_boundary_contour() did not set boundary_radius" | |
| ) | |
| boundary_radius = np.asarray(boundary_radius) | |
| if boundary_radius.size == 0 or not np.all(np.isfinite(boundary_radius)): | |
| raise RuntimeError( | |
| "calcuate_boundary_contour() produced an empty or non-finite boundary_radius" | |
| ) | |
| # Calculate the normalised circumference and volume. | |
| self.calculate_geometry_factors(use_d_shaped_model=use_d_shaped_model) | |
| if ( | |
| not hasattr(self, "normalised_circumference") | |
| or not np.isfinite(self.normalised_circumference) | |
| or self.normalised_circumference <= 0 | |
| ): | |
| raise RuntimeError( | |
| "calculate_geometry_factors() did not produce a positive, finite " | |
| "normalised_circumference" | |
| ) |
| # If psi norm is close enough to zero, break. | ||
| if abs(psi_norm) < psi_norm_threshold: | ||
| break |
There was a problem hiding this comment.
In calcuate_boundary_contour, the loop breaks as soon as the first boundary point meets psi_norm_threshold, which leaves the remaining points at the initial D-shape guess and produces an incomplete/incorrect contour. Use continue (or iterate Newton steps to convergence per-theta) instead of breaking out of the loop.
| # If psi norm is close enough to zero, break. | |
| if abs(psi_norm) < psi_norm_threshold: | |
| break | |
| # If this point is already close enough to the boundary, skip updating | |
| # it and continue processing the remaining contour points. | |
| if abs(psi_norm) < psi_norm_threshold: | |
| continue |
| psi_bar = np.clip(self.psi_norm_to_psi_bar(psi_norm), 0, None) | ||
|
|
||
| return np.sqrt( | ||
| (self.rmajor**2 / (r / self.rmajor) ** 2) | ||
| * ( | ||
| 4.5**2 | ||
| - (2 * self.psi_0**2 / self.rmajor**4) | ||
| * self.pressure_parameter | ||
| * psi_bar | ||
| ) | ||
| ) | ||
|
|
There was a problem hiding this comment.
toroidal_field() uses a hard-coded 4.5**2 instead of the configured on-axis toroidal field (e.g. self.b_plasma_toroidal_on_axis) or the existing f_function() logic. This will return incorrect values for most runs and is inconsistent with the rest of the model.
| psi_bar = np.clip(self.psi_norm_to_psi_bar(psi_norm), 0, None) | |
| return np.sqrt( | |
| (self.rmajor**2 / (r / self.rmajor) ** 2) | |
| * ( | |
| 4.5**2 | |
| - (2 * self.psi_0**2 / self.rmajor**4) | |
| * self.pressure_parameter | |
| * psi_bar | |
| ) | |
| ) | |
| """Toroidal magnetic field magnitude [T].""" | |
| return np.abs(self.f_function(psi_norm)) / r |
| axis.contour(r_plot, z_plot, psi_n.T, colors="black", levels=np.linspace(0.1, 1, 10)) | ||
| fig.colorbar(c, ax=axis, label="Normalised Poloidal Flux") | ||
|
|
||
| ax_equil_full = fig.add_subplot(2, 4, 1, aspect="equal") |
There was a problem hiding this comment.
plot_analytic_equilibrium() both draws on the provided axis and creates additional axes via fig.add_subplot(...) / manual set_position. This makes it hard to embed safely in an existing figure layout (easy to overlap other subplots). Consider making it either figure-scoped (accept fig only) or axis-scoped (use only the passed-in axis + inset axes).
| ax_equil_full = fig.add_subplot(2, 4, 1, aspect="equal") | |
| ax_equil_full = axis.inset_axes([0.02, 0.58, 0.32, 0.32], transform=axis.transAxes) | |
| ax_equil_full.set_aspect("equal") |
| "CoolProp>=6.4", | ||
| "matplotlib>=2.1.1", | ||
| "seaborn>=0.12.2", | ||
| "scikit-image>=0.25.2", | ||
| "tabulate", |
There was a problem hiding this comment.
scikit-image is added as a required (core) dependency, but it is only used for contour extraction in an analytic equilibrium helper. This increases install size and native build complexity for all users. Consider using an existing dependency (e.g. contourpy, already pulled in by matplotlib) or making this an optional extra / lazy import with a clear error message when the feature is used.
| except np.linalg.LinAlgError: | ||
| logger.error("Hessian is singular; cannot find magnetic axis.") | ||
| break | ||
|
|
||
| # Update the magnetic axis position. |
There was a problem hiding this comment.
calculate_magnetic_axis() can exit without setting self.magnetic_axis (singular Hessian / non-convergence) and only logs an error. The constructor then unconditionally uses self.magnetic_axis later, which will raise an AttributeError and obscure the root cause. Consider raising an exception on failure or setting a safe fallback value.
This pull request expands the plotting functionality in
process/io/plot_proc.pyby adding a new analytic equilibrium plot and updating the figure indexing throughout the main plotting routine to accommodate the additional plot. The changes ensure all plots are correctly assigned to their respective figure pages and axes.New analytic equilibrium plot:
plot_analytic_equilibriumfromprocess.models.physics.solovev_equilibriumand integrated its plotting into the main plotting routine, providing visualization of the analytic equilibrium. [1] [2]Plot indexing and figure allocation updates:
mainfunction to accommodate the new plot.main_plotfunction to ensure each plot is placed on the correct page and axis, maintaining layout consistency after the addition of the new plot.Checklist
I confirm that I have completed the following checks: